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1  Introduction 

Multitarget  Multisensor  Tracking  (MTMST)  as  the  name  implies,  is  the  process  of 
tracking  a  multiple  of  targets  using  a  number  of  sensors.  In  general  terms  this  process  is 
divided  into  two  major  sub  processes,  one  has  to  do  with  the  data  estimates  at  each  discrete 
point  in  time,  and  the  second  with  the  development  of  the  target  trajectory  based  on  the 
data  estimates.  The  process  of  combining  data  from  a  number  of  sensors  taken  at  a  discrete 
point  in  time  is  called  sensor  fusion,  and  the  corresponding  combinatorial  problem  is  the 
data  association  problem. 
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When  tracking  a  single  target  with  a  single  sensor,  noisy  measurements  arc  made  at 
discrete,  periodic  points  in  time  from  the  sensor  regarding  the  moving  target.  The  goal  is  to 
establish  a  time- continuous  estimate  of  the  target’s  trajectory.  For  the  MTMST  problem, 
noisy  measurements  arc  made  from  an  arbitrary  number  of  spatially  diverse  sensors  regard¬ 
ing  an  arbitrary  number  of  targets  with  the  goal  of  estimating  the  trajectories  of  all  of  the 
targets  present.  Furthermore,  the  number  of  targets  may  change  (e.g.,  by  wandering  out 
of  and  into  detection  range)  during  the  problem  time  interval  and  the  sensors  themselves 
may  be  in  motion.  For  example  in  Figure  1,  two  sensors  observe  three  targets,  thereby 
wc  have  three  measurements  from  each  sensor.  The  MTMST  problem  is  to  partition  all 
six  measurements  (three  from  each  of  the  2  sensors)  into  three  subsets  that  correspond  to 
targets. 


rag  replacements 

x 

V 

sensor  1 
sensor  2 
terrain 
targets 

Figure  1:  Sensor  fusion 

Target  state  estimates  for  single-target  single-sensor  tracking  arc  obtained  satisfactorily 
by  applying  one  of  several  variants  of  the  Kalman  filter.  However,  no  method  known  can 
consistently  obtain  good  trajectory  estimates  for  the  more  general  MTMST  problem. 

The  MTMST  has  generally  been  divided  into  2  sub  problems:  state  estimation/propagation 
and  data  association.  The  combinatorial  nature  of  the  MTMST  problem  results  from  the 
data  association  sub  problem;  that  is,  given  that  each  sensor  m,  m  =  1 . . .  M ,  makes  nm 
measurements,  n*  /  how  do  we  correctly  partition  the  measurements  into  all  of 

the  targets,  ensuring  that  all  of  the  measurements  arc  used  exactly  once?  Total  enumera¬ 
tion  of  all  possible  partitions  is  obviously  unrealistic,  since  the  size  of  the  search  space  grows 
exponentially. 

In  MTMST  a  data  association  problem  is  presented  at  discrete,  periodic  moments  in 
time.  Wc  expect  that  by  obtaining  the  correct  partition  or  association,  the  data  from  all 
of  the  sensors  may  be  combined  in  some  manner  to  make  estimates  of  target  positions  that 
arc  more  accurate  than  that  obtained  from  a  single  measurement. 

As  it  is  developed  in  the  paper  by  Murphcy,  Pardalos  and  Pitsoulis  [9] ,  the  data  associ¬ 
ation  problem  with  respect  to  the  MTMST  can  be  modeled  mathematically  as  a  Multidi¬ 
mensional  Assignment  Problem  (MAP).  In  this  research  wc  will  examine  the  following 
aspects  of  the  MAP: 

•  Problem  Generators  that  will  produce  problem  instances  with  known  optimal  solu¬ 
tions. 

•  Efficient  data  structures  to  handle  the  massive  data  sets  associated  with  the  instances 
of  the  MAP. 
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•  Heuristic  algorithms  for  providing  near  optimal  solutions  to  the  MAP. 

•  Lower  bounds  and  exact  algorithms  for  solving  the  MAP. 

•  A  comprehensive  survey  on  the  MAP  which  will  report  all  the  latest  research  devel¬ 
opments  for  this  important  combinatorial  optimization  problem. 

The  objective  of  this  research  is  to  develop  efficient  algorithms  for  solving  the  multidi¬ 
mensional  assignment  problem  in  real  time. 


2  Formulations 

For  most  combinatorial  optimization  problems,  different  mathematical  formulations  exist 
that  are  equivalent.  They  usually  reflect  the  effort  to  view  the  structural  characteristics  of 
the  problem  from  a  different  angle,  in  hope  of  developing  new  solution  techniques.  In  this 
section  we  will  present  the  various  equivalent  mathematical  formulations  for  the  MAP.  Let 
us  denote  a  multidimensional  assignment  problem  of  dimension  m  with  MAP(m). 

2.1  0-1  Integer  Programming  Formulation 

In  what  follows  the  optimization  problem  formulation  of  the  multidimensional  assignment 
problem  will  be  presented,  first,  as  an  Integer  Program  (IP).  Let  ni5  n2, . . . ,  nrn  be  integers 
which  define  the  ranges  in  each  dimension.  The  input  data  consists  of  one  cost  array 
C  e  Mni  XTl2X  ‘xn'mJ  where  without  loss  of  generality  assume  that  n\  <  n-2  <  •  *  *  <  nm.  The 
IP  formulation  of  the  m-dimensional  assignment  problem  is  the  following: 

Y,  ch...imxh...im  (i) 

»1=1  im  —  1 

n-2  nin 

y  ^  ^  V  "‘im  =  lj  U  1 }  *  *  ’  3  ^1  )  (^) 

12  =  1  im=l 

ni  rik~  1  nk+1  nm 

X  "  X  X/  ■  '  '  X/  Xh-im  ^  (3) 

tl  =1  ik -1=1  ifc+i==l  tm= 1 

for  ik  =  1, . . . ,  rik  and  k  =  2, . . . ,  m  -  1  (4) 

ni  Tlm-l 

y  *  *  ’  y  ]  %ii  •  im  —  Ij  ^rn  =  1?  *  •  •  i  (^) 

*1=1  i,n-l=  1 

®ii  ...im  €  {0, 1}  for  all  ii  •  •  •  im.  (6) 

We  will  call  the  constraints  in  the  above  formulation  multidimensional  assignment  con¬ 
straints.  The  inequalities  arc  a  direct  consequence  of  the  fact  that  the  ranges  ,  n2, . . . ,  nm 
arc  not  equal.  However  we  could  assume  without  loss  of  generality  that  m  =  ri2  =  •  •  •  =  nm, 
by  introducing  artificial  costs  in  the  cost  array  C  €  KniXn2X",xnm  with  a  very  high  value 
(thereby  forcing  them  not  to  be  included  in  the  optimal  solution) . 

Example  2.1  Consider  that  we  have  m  =  3,  or  in  other  words  we  have  the  well  known 
three  dimensional  assignment  problem  where  using  the  0-1  formulation  stated  in  (2)- (6)  we 


MAP/p (m)  :  min 

s.t. 
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will  have: 

711  Tl2  Tl3 

MAPjp(3)  :  min  £]EE  CijkXijk 

i= 1  j= 1  fc=l 

712  713 

EE  X^ijk  —  1)  ^  —  !)•••?  ^1  j 

j=l  k=  1 

71 1  713 

^  ^  &ij  J  !>•••>  ^2 1 

i=l  Jt=l 

Til  712 

EE  %ij  k  ^  1  j  k  —  !»•••?  ^3  ? 

i=l  j‘=l 
%ij 

where  in  the  above  formulation  we  assume  that  n i  <  ri2  <  713. 

The  above  mentioned  IP  formulation  of  the  MAP  however  is  not  very  convenient  with  respect 
to  algorithmic  issues,  it  does  however  provide  an  insight  to  the  nature  of  the  problem.  For 
example  it  can  be  used  to  derive  lower  bound  for  the  problem  by  relaxing  the  integrality 
constraints  and  solving  the  linear  program.  It  has  also  been  used  in  a  Lagrangian  relaxation 
solution  scheme  in  [13]. 


2.2  Permutations 

In  this  section  we  will  offer  a  different  formulation  for  the  MAP  using  permutations.  This 
formulation  will  provide  essential  insight  so  as  to  construct  efficient  data  structures  for  the 
problem  and  local  search  procedures.  The  main  characteristics  of  this  formulation  is  its 
compactness,  and  its  combinatorial  nature. 

Let  Xm  be  the  set  of  all  arrays  X  =  (®i1...im)n1x-xnm  such  that  their  elements  satisfy 
the  constraints  (2)  through  (6).  We  will  call  Xm  permutation  arrays. 

Given  a  set  of  integers  N  =  {1,  2, . . . ,  n}  a  partial  ^-permutation  on  N  will  be  defined  as 

/  1  2  •  •  •  k  \ 

the  mapping  p  :  {1,2, .. . ,  fc}  N,  where  k  <  n.  We  write  p  :=  [  ^  ^  ^  J 

but  for  short  from  now  we  will  write  p  =  (p(l),p(2), . . .  ,p(fc)).  For  example  for  N  = 
{1,2, 3, 4, 5}  a  partial  3-permutation  might  be  p  =  (5,3,1).  Denote  the  set  of  all  fc- 
permutations  of  a  set  of  integers  N  =  {1, 2, . . .  ,n}  as  n^fc)- 


Remark  2.1  Given  that  |AT|  =  n  and  some  integer  k  =  1,2,  ...,n  we  have  |II(Ar,/c)|  = 

n\ 

(n  — fc)!  * 

Define  the  sets  of  integers  Ni  :={1,2 i  =  1, 2, . . . ,  m,  and  let  be  the  set 

of  all  ni -permutations  pi  :  N\  — ►  N{.  Note  that  |ri( iv£,ni) I  =  (ni  ~n7)i  >  f°r  eac^  *  =  1, . . .  ,m. 
Given  pi, . . .  ,pm  let  us  define  a  matrix  P  :=  (pi, . . .  ,pm),  that  is  P  can  be  regarded  as  a 
n\  x  m  matrix  with  columns  the  n\ -permutations  pi, . . .  ,pm, 


/Pl(l)  PiiX)  •••  Pm(l)\ 

Pi  (2)  P2(2)  •••  Pm  (2) 

\Pl(«l)  P2(m)  •••  Pm(ni)/ 
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Now  given  n\  <  n2  <  *  •  •  <  nm,  and  denoting  by  =  (1, 2, . . . ,  nx)  the  identity  permuta¬ 
tion  of  N\ ,  let  Pm  be  the  following  set 

Pm  :=  { P  =  (pi  ;  Pi  =  ejVi ,  Pi  €  nWini),i  =  2, ...  ,m}. 

Remark  2.2  |P-|  =  1^ 

Note  that  there  is  a  one-to-one  correspondence  between  the  elements  of  Xm  and  Pm  for 
given  ni  <  ri2  <  ■  ■  ■  <  nm,  since  given  any  P  G  Pm  we  can  construct  an  X  G  Xm  as 

_  fl  if  Pj (^i)  =  V? .?  =  2, . . . , ?n5  /q\ 

11  lm  1^0  otherwise, 

for  ii  =  1, . . .  ,m,  and  vice  versa.  We  can  therefore  formulate  the  MAP  using  permutations 
as  follows: 


MAPn(m)  : 

min 

n  1 

^  ^iP2(0*“ Pm-l(0Pm(0 
i=l 

(9) 

s.t. 

Pi  ^  TE(/\/'i  ,n.x )  ’  ^  2,  .  . 

.  ,?7l. 

A  solution  to  the  MAP  using  the  above  formulation  will  be  any  P  G  Pm  such  that  the  above 
objective  function  is  minimized.  Moreover  from  remark  2.2  we  can  conclude  about  the  size 
of  the  feasible  space  which  is  finite,  but  of  course  it  grows  factorially  with  the  problem  size 
and  dimension. 

Example  2.2  Consider  that  we  have  m  =  2  with  |Ai|  =  ni  and  |7V2|  =  n2  with  ri\  <  n2. 
This  is  the  well  known  linear  assignment  problem : 

m 

min  Cip(0  (10) 

i=  1 

S.t.  p  G  n(iV2,ni)» 

where  it  is  well  known  that  it  can  be  solved  in  0(n3)  computational  time. 

2.3  Graph  Theoretic 

Consider  a  three  dimensional  assignment  problem  where  we  have  G  =  (V,  E)  and  V  = 
V\  U  V2  U  V3.  The  graph  will  contain  edges  and  hyperedges ,  and  the  edge  set  will  be  composed 
as  follows,  E  =  E\  U  P2  U  £3  U  E4  where  E\  =  V\  x  V2,  E2  =  V2  x  V3,  E$  =  V\  x  V3  and 
finally  E4  =  V\  x  V2  x  V3.  We  will  denote  an  element  of  P  by  e#  if  it  is  the  edge  joining 
vertices  i  and  j,  or  eijk  if  it  is  the  hyperedge  joining  the  vertices  i,j  and  k.  An  instance  of 
such  a  three-partite  graph  with  a  feasible  solution  is  given  in  Figure  2.  Based  on  the  above 
example  we  can  now  formulate  the  multidimensional  assignment  problem  as  the  following 
problem  on  graphs: 

MAPc?(m):  Given  an  m-partite  graph  G  —  (V,  E)  with  vertex  set:  V  =  V\  U  V2  U  •  ■  •  U  Vrn 
were  | Vi|  <  \V2\  <  *  •  *  <  \Vrn\,  and  a  set  of  hyperedges  E  C  V\  x  V2  x  •  •  •  x  Vm,  where  equality 
holds  when  the  problem  is  dense ,  and  a  cost  function  c  :  E  — ►  R,  find  \  Vi\  disjoint  cliques 
of  size  m  with  minimum  sum  of  edge  costs. 
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PSfrag  replacements 
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Figure  2:  A  tri-partite  graph 

2.4  Inner  Product 

This  formulation  of  the  MAP  is  an  immediate  result  from  the  integer  programming  formu¬ 
lation  of  the  problem,  by  making  the  following  generalization  of  the  inner  product.  More 
specifically,  let  the  inner  product  between  matrices  A,  B  €  MmXn  be  defined  as 

m  n 

( A ,  B)  aij bij • 

i= 1  j=l 

We  can  generalize  the  above  for  higher  dimensional  arrays.  Define  the  inner  product  between 
two  arrays  A,B  e  to  be 

n\  nm 

(A,  B)  :=  a»l— 

i1=l  zm  — 1 

Therefore  we  can  now  formulate  the  MAP  as: 

MAPo(m)  :  min  (C,  X)  (11) 

s.t.  x  e  xm, 

where  C  is  the  m-dimcnsional  cost  array,  and  Xm  the  set  of  all  m-dimcnsional  permutation 
arrays. 

3  Matroid  Intersection 

The  multidimensional  assignment  problem  (MAP)  can  be  formulated  as  the  intersection  of 
partition  matroids  as  it  will  be  shown  in  this  section.  A  whole  section  is  devoted  to  the  sub¬ 
ject,  since  this  formulation  is  more  general  than  all  the  previously  mentioned  formulations, 
and  it  presents  a  new  way  of  looking  into  the  problem  which  introduces  many  interesting 
open  questions  regarding  the  problem. 
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3.1  Definitions 

Lot  us  first,  state  some  necessary  definitions  from  matroid  theory  [14]. 

Definition  3.1  Given  a  finite  set  E  and  some  family  of  subsets  F  C  2E  ,  the  set  system 
M  =  ( E ,  F)  is  called  a  matroid  if  the  following  axioms  hold: 

(II.)  %eF. 

(12.)  If  X  £  F  and  Y  C  X  =*  V  G  F. 

(13.)  IfX,YeF  and  \X\  =  \Y\  +  1  =►  3  x  €  X/Y  such  that  Y  U  {x}  e  T. 

If  only  (II)  and  (12)  are  true  then  the  system  is  called  independent  system. 

Axiom  (12)  will  be  called  the  subinclussiveness  property,  while  axiom  (13)  the  accessibility 
property.  The  elements  of  IF  arc  called  independent,  while  those  not  in  T  arc  called 
dependent.  Maximal  independent  sets  are  called  bases  and  minimal  dependent  sets  are 
called  circuits.  Given  some  independent  system  ( E ,  T ),  the  rank  of  some  /  C  E  is  defined 
as 

r(I)  :=  max{|A|  :  X  C  I,X  £  F}1 

while  the  lower  rank  as 

lr(I)  :=  min{|X|  :  X  C  7,  X  €  F,  X  U  {s}  $  F  Vs  €  I/X}. 

We  can  interpret  the  values  of  r(E)  and  lr(E)  as  the  cardinality  of  the  maximum  cardinality 
basis  and  the  minimum  cardinality  basis  of  the  independent  system  respectively.  Note  that 
in  the  case  where  the  independent  system  is  a  matroid,  then  r(E)  =  lr(E).  For  any  matroid 
M  =  ( E ,  F)  there  is  an  associated  independent  system  on  the  same  ground  set  M*  =  (F,  F*) 
called  the  dual  of  M,  such  that  the  bases  of  M*  are  the  complements  of  the  bases  of  M. 
For  independent  systems  in  general,  the  dual  of  an  independence  system  (E,F)  is  defined 
as (E,F*)  where 

F*  :=  {/  C  E  :  there  is  a  basis  B  of  (E,  F)  such  that  I  n  B  =  0} 

3.2  Formulation  of  the  MAP(m) 

In  order  to  state  the  matroid  intersection  formulation,  we  will  construct  a  hypergraph  asso¬ 
ciated  with  the  problem.  Given  ni,ri2, . . .  ,nm  as  above  construct  the  m-partite  hypergraph 
H(Vy  E)  as  follows: 

m 

•  V  :=  |J  Vi ,  where  |V5|  =  n*  and  Vi  n  Vj  =  0  for  all  i,  j. 

i~l 

•  E  :=  {(v<1}  Vi2, . . .  yVim)  :  Vi.  e  VjJ  =  1, . . .  ,m},  or  in  other  words,  E  consists  of  only 
hyperedges  on  V  of  “size”  m  only.  Equivalently  we  can  define  E  :=  V\  x  •  •  •  x  Vm. 

•  A  cost  function  c  :  E  —*•  R,  which  is  defined  by  the  cost  array  C.  That  is,  given  any 

hyperedge  e  =  (u^ ,  , . . .  ,u*m)  the  associated  cost  is  c(e)  = 

Given  any  two  hyperedges  e  —  (v^ ,  Vi2 , . . . ,  virn )  and  w  =  (Uj: ,  Vj2 , . . . ,  Vjm )  we  say  that  they 
arc  disjoint  iff  Vit  ^  Vjt ,  /  =  1, . . . ,  m.  In  the  MAP(m)  we  have  to  find  ni  disjoint  hyperedges 
of  minimum  total  weight.  In  other  words  we  have  to  find  a  maximal  subset  I  C  E  of  disjoint 
hyperedges,  with  minimum  weight  c(I)  :=  Y)  c(e)-  If  problem  is  dense,  which  implies 
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that  the  multidimensional  array  C  is  dense,  then  the  cardinality  of  the  maximal  subsets  I 
would  be  equal  to  n\ . 

Now  consider  the  following  partitions  of  E  as  defined  by  Vi, . . . ,  Vm.  Let 

UVj  =  u  PM,v)j  j  = 

veVj 

be  m  partitions  of  E ,  where  for  each  v  £  Vj  we  have  the  following  set  of  hyperedges  defined 
by  it: 

P(Vj,v)  {e  €  E  :  e  is  incident  to  vertex  v  £  Vj }. 

When  we  say  that  a  hyperedge  e  =  (v^ , . . . ,  V{m )  is  incident  to  some  v  £  Vj  we  mean  that 
V(j  =  v.  If  we  define  the  independence  relationship  as  the  following  family  of  subsets: 

TVj  :=  {/  C  E  :  \lnP{VjyV)\  <  l,for  all  v  £  Vj},  j  =  1,.. .  ,m. 

In  words,  Ty3  is  a  family  of  subsets  of  the  hyperedges  of  H  such  that  no  subset  contains 
any  two  hyperedges  which  arc  incident  to  the  same  vertex  among  those  of  Vj. 

Theorem  3.1  Given  an  MAP(m),  the  set  system  Myj  =  (E,Tyj)  is  a  matroid,  j  = 

Proof  3.1  Axioms  (II )  and  (12)  are  evidently  true.  For  (13),  assuming  that  we  have  X,  Y  £ 
Tv.  with  \X\  —  \Y\  +  1,  we  have  to  show  that  3e  €  X/Y  such  that  Y  U  {e}  €  Ty5 .  Note 
that  the  edges  of  any  given  set  I  C  E,  due  to  the  independence  relationship,  determine  \I\ 
distinct  sets  P(yj,v)-  V  any  two  were  not  distinct  that  would  imply  that  there  exists  some 
vertex  v  £  Vj  such  that  \I  n  P(vjjV)\  =  2.  Therefore  since  \X\  =  \Y\  +  1,  there  exists  some 
vertex  v  G  Vj  such  that  \X  D  P(V}.  v)|  =  1  and  | Y nP(yj?v)|  =  0,  or  in  other  words  there  exists 
an  edge  e  £  X/Y  such  that  Y  U  {e}  £  Ty. .  □ 

In  particular  the  independence  system  above  is  a  partition  matroid.  In  the  MAP(m) 

m 

we  want  to  find  some  I  £  f]  Ty.  of  maximum  cardinality  and  minimum  weight  c(/).  This 

i=i 

is  a  weighted  matroid  intersection  problem  which  will  we  define  as  follows: 


MAP^(m):  Given  m  matroids  My.  =  (E,Ty.)  by  their  independent  oracles, 

m 

and  a  cost  function  c  :  E  — >  M,  find,  a  set  I  £  p)  Ty5  of  maximum  cardinality 

i=j 

such  that  its  weight  ^2  c(e)  its  minimum, 
eei 


rn 

Note  that  the  set  system  ( E ,  p|  Ty5)  is  not  a  matroid,  the  intersection  of  matroids  is  not 

j=i 

matroid  in  general  (it  is  easy  to  give  a  counterexample  for  the  specific  system  that  shows 
that  its  not  a  matroid).  We  do  have  however 

m 

Corollary  3.1  ( E ,  f]  Ty 7.)  is  a  independence  system. 

3=1 

It  follows  directly  using  theorem  3.1,  since  the  intersection  of  matroids  is  an  independence 
system  [14]. 
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algorithm  independent  (/) 

1  Ve  =  {vi1,Vi2,...,vim)  €  /  do  -*■ 

2  do  j  =  1, 2, . . .  ,m  — ► 

3  if  is  labeled  STOP  (7  £  ^.) 

4  else  label  ; 

5  od; 

6  od; 

7  return  7  €  ; 

end  independent; 


Figure  3:  Independence  Oracle 


3.3  Oracles 

A  word  has  to  be  said  here  regarding  the  independence  oracles  mentioned  in  the  problem 
definition.  An  independent  oracle  can  be  considered  as  an  algorithm  that  takes  as  an  input  a 
set  7  C  and  answers  the  question  whether  this  set  is  independent  or  not  (i.e.  that  is  satisfies 
the  independence  relationship).  In  our  case  it  is  trivial  to  see  that  given  a  set  of  hyperedges 
I  C  E  and  some  j  e  {1, . . . ,  m},  we  can  check  whether  7  €  Tyj  in  0(\I\)  time  by  examining 
whether  there  arc  any  two  edges  in  7  with  their  jth  component  the  same.  Similarly  we  can 

rn 

check  in  linear  time  whether  I  e  f]  TVi .  The  following  simple  labeling  algorithm  shown  in 

i=i 

Figure  3  can  be  considered  as  an  independent  oracle. 

It  is  clear  that  the  algorithm  will  require  0{k\I\)  computational  time.  Checking  indepen¬ 
dence  of  a  set  is  needed  when  the  algorithm  used  for  constructing  the  solution  proceeds  in 
an  additive  fashion,  starting  from  the  empty  set  and  adding  one  element  at  a  time  until  the 
set  is  independent.  If  we  wish  to  proceed  in  the  other  way,  that  is  starting  from  the  ground 
set  E  and  subtracting  elements  from  it  until  we  reach  a  solution,  we  need  a  different  type 
of  oracle,  a  basis  oracle.  A  basis  oracle  can  be  considered  as  an  algorithm  which  takes  as  an 
input  a  set  7  C  E  and  gives  an  answer  to  the  question  of  whether  or  not  this  set  contains  a 
basis.  The  construction  of  such  a  basis  oracle  algorithm  remains  an  open  question. 


3.4  Bounds 

m 

Establishing  the  fact  that  the  system  ( E ,  p|  J~Vj)  is  independent  enables  us  to  use  existing 

j- 1 

results  on  bounds  for  greedy  type  algorithms  for  these  type  of  problems.  Specifically,  given 
any  (E,T)  be  an  independence  system  and  some  cost  function  c  :  E  — ►  R+,  a  greedy  type 
algorithm  will  produce  a  solution  to  the  minimization  problem  by  starting  with  the  the 
ground  set  E  and  start  excluding  elements,  always  taking  out  the  element  with  the  largest 
weight.  The  algorithm  will  stop  when  the  subset  at  hand  is  a  basis,  a  fact  which  can  be 
checked  using  the  basis  oracle  mentioned  in  the  previous  section.  The  following  result  is  due 
to  [4]. 

Theorem  3.2  (Korte  and  Monma  [1979])  Let  (E,E)  be  an  independence  system .  For 
c  :  E  — ►  M+  let  G(E,E,c)  denote  a  solution  found  by  a  greedy  type  algorithm ,  while 
OPT{E ,  c )  is  the  optimum  solution .  Then 


1  <  G(E,E,c) 

~  OPT{E,F,c) 


< 


\F\-lr*(F) 

max  T— - T-=rr- 

fce  |F|  -  r*(F) 


p{E,T) 
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for  all  c  :  E  — >  M+J  tu/iere  r*(F)  and  lr*(F)  are  the  lower  rank  and  rank  respectively  as 
defined  on  the  dual  (E^J7*).  There  is  a  cost  function  where  the  upper  bound  is  attained. 

Based  on  the  above  theorem,  we  can  see  that  if  we  can  evaluate  the  value  of  p(E,  J7)  we  can 
establish  a  theoretical  bound  on  the  performance  of  a  greedy  type  algorithm.  Note  however 
that  the  above  result  concerns  all  independent  systems,  and  it  is  an  open  question  whether  or 
not  for  the  independence  system  established  for  the  multidimensional  assignment  problem, 
there  is  a  closed  form  expression  for  p(E,  J7),  which  in  turn  would  imply  an  ^-approximation 

771 

algorithm  for  the  problem.  It  is  conjectured  that  we  have  p(E,  fj  J7yj )  =  /(ra),  that  is  the 

3= i 

value  of  p  is  some  constant  function  on  the  dimension  of  the  problem  m. 


4  Data  Structures 

It,  immediately  follows  from  the  nature  of  the  problem  that  an  efficient  way  of  handling 
the  large  data  is  required.  There  arc  two  basic  problems  in  handling  the  input  data  of  an 
instance  of  the  MAP.  First,  multidimensional  arrays  arc  not  desirable  in  terms  of  memory  re¬ 
quirements,  and  second,  the  variable  dimension  m  poses  implementation  problems.  Both  of 
these  problems  are  solved  if  we  transform  the  multidimensional  array  into  a  one-dimensional 
array. 


4.1  Transformation  into  one  dimension 

Let  hi  nLin*  for  *  =  1 ,  The  array  C  G  riXn2x-xnm  js  stored  as  an  one- 

dimensional  array  c  G  R*,  in  column  format.  For  example  if  C  G  R2x2x2  then  the  array 
c€l8  will  be 

c  =  (Cm,  C21U  C121,  C22U  C 112,  C2125  £122?  C222)' 

In  addition  to  the  sets  of  integers  TV*,  i  =  1, . . . ,  m  defined  previously  define  also  the  set 
N  :=  {1, 2, . . . ,  7r}.  The  set  N\  x  ■  ■  ■  x  Nm  can  be  regarded  as  the  index  set  of  C  while  N 
the  index  set  of  c.  The  transformation  between  an  index  {i\,i 2, . . . ,  im )  of  C  to  an  index  k 
of  c  is  given  by  the  following  function  /  :  Ni  x  •  *  -  x  Nm  — +  N 

771 

f  (H  >  ^2l  •  *  *  5  i 771 )  * —  ^1  T  'y  ^  "Rj  —  l  (fj  l)  (^) 

3=  2 


The  inverse  (from  an  index  k  of  c  to  an  index  (fi,Z2, . . .  ,  fm)  of  C)  can  be  obtained  in  a 
similar  fashion  from  the  function  f~l  :  N  — ►  Nx  x  ■  •  •  x  Nrn,  the  inverse  of  /,  which  is 
defined  by 


}~\k)  := 


(  |fc  -  1  -  YT=l  'Krn-l{im-l+l  ~  1)J  +  1  \ 

I  k~l~T7Llq  +  I  ,  y 

[_  7Tm  —  1  J 

[SJ+1 


(  H  \ 


V  im  J 


Lemma  4.1  The  function  f  is  a  bijection. 
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Proof  4.1  Since  \  N\  x  •  •  •  x  JVm|  =  | N\  =  7r  by  construction ,  for  f  to  be  bijective  it  suffices  to 
show  that  for  any  two  (ii, . . .  ,im),  (ji, . . . ,  jm)  €  Ni  x  •  •  ■  x  A^m,  (ii, . . .  ,im)  7^  (ji»  •  ■  • >Jm) 

(i.e.  distinct  elements),  we  have  /(ii,...,im)  /  /(ju  •  •  •  ,  jm)-  Clearly  this  is  true  for 

m  -  l,  since  f(i{)  =  ii  ^  ji  =  /(ji).  Suppose  that  is  true  for  m  and  let  (ii, . . . , im  +  1)  ^ 
(ji,  •  •  • >  jm.  +  1)-  Then  we  have 

f{i  1  j  *  •  •  j  im  +  1)  —  /(il ,  .  .  ■  }  ^m)  H“  TTm (im+1  ~  1) , 

f{jli  •  *  •  j  Jm  +  1)  =  /  ( J 1  j  •  ■  ■  ?  Jm)  H"  “  !)• 

T/tere  are  three  cases:  i)  7^  (ji>  •  •  •  1  jm)  and  im+i  =  jm+i>  (ii,...,im)  = 

(ii*  *  •  * >  Jm)  and  therefore  im+i  ^  jm+ 1,  and 

Hi)  (ii, . . .  ,tm)  7^  (ii*  *  ■  •  1  Jm)  and  im+i  ^  jm+ 1-  /n  a//  0/  i/ie  cases  /(ii, . . .  ,im  +  1)  7^ 
/(ii, . .. ,  jrn  +  1)  since  we  h,a.ve  respectively  for  the  first  case  the  induction  hypothesis  a,nd 
for  the  second  case  im+i  #  jm+i*  Tor  the  third  case,  note  that  f(ii, . . .  ,im),  / (ii , . .  - ,  Jm)  € 
[1, . . . , 7rm]  and  7rm(im+i  -  1),  7rm(jm+i  -  1)  €  [0, 7rmj . . . ,  (nm+i  -  l)ir m],  which  imply  the 
following  inequalities  considering  that  we  have  distinct  elements 

1  ^  |/(il  ?  *  •  •  j  ^m)  “  / (jl  i  *  •  ■  i  Jm)  I  ^  ~  h 

[7T.m(im+l  “  1)  ^m(jm+l  ”  1)1  —  ^m* 

Therefore,  we  have  that 


‘7Tm(im+ 1  l)  ^  7Tm(jm+x  l)  +  (/(jl ,  •  •  •  ,  Jm)  /(il ,  . .  .  ,  im))j 


wtoic/i  implies  that  f{i\ , . . . ,  im  +  1)  ^  /(ji, . . . ,  jm  +  1). 


□ 


From  /  and  /  1 


it  is  clear  that  wc  can  construct  the  c  array  as  c(i)  =  C(f  1(i)),i  = 


4.2  The  index  structure  of  c 

Assume  that  wc  arc  given  an  element  (ii,  i2, . . . ,  im)  of  C,  and  for  a  given  fc  we  wish  to 
generate  the  set 


iV  3  Sik  i —  {i  £  iV  :  /  (i)  —  (ii , . . . ,  , . . . ,  im ) }  5 

that  is  the  set  of  all  elements  in  N  such  that  their  corresponding  element  in  Nx  x  •  •  •  x  Nrn 
contains  ifc.  A  direct  approach  would  be  to  generate  the  set  Sik  C  N\  x  •  ■  *  x  A^m  defined  as 

S^k  { (ii  ?*•  •  jifc?  ’  •  *  5  im)  *  ij  —  h  *  *  *  ?  ,  j  1} . . . ,  m,  j  7^  , 

and  then  take  the  image  of  S'h  under  /,  Sik  =  f(S'ik ),  a  process  which  will  require  D(m ^) 
time.  The  way  that  c  is  constructed,  it  follows  that  its  index  structure,  that  is  the 
relationship  between  the  indices  of  c  and  those  of  C ,  is  well  defined.  Let  us  define  for 
expository  purposes,  an  ordered  7r-tuplc  which  corresponds  to  the  index  set  N  of  c,  as 
I  :=  (1, 2, . . . ,  7r  -  1, 7r).  For  every  i  e  I  there  is  a  corresponding  element  (ii,  i2, . . . ,  im)  in 
the  index  set  of  C.  Define  the  following  for  every  dimension  k  =  1,  2, . . .  ,m  and  element 
ik  6  Nk : 

*=  (ife  1  )tTA:  —  1  lj 
7 r 
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According  to  the  above  values,  for  any  4  we  can  identify  the  following  irk-t -sub-tuples 
contained  in  J, 


/  =  (1,2, . . .  ,B\k, . .  . .,  £**, . 1,tt), 

7Tfc  (l^fc“3)7rfc 


where  =  ((Z-l)7Tfc  +  0ifcI  (/-  l)7Tfc  +  aifc  +  1, . . . ,  (/-  l)7rk  +  Vik  +7Tfc_i  - 1)  for  l  —  1, . . . ,  v^. 
So  for  every  4»  I  contains  sub-tuples  of  size  7Tfc-i,  starting  from  the  clement  aik  and 
occurring  every  nk  elements. 


Lemma  4.2  Sik  =  IJ^i  • 

Proof  4.2  5mce  £4  =  f(Sik)  it  suffices  to  show  that  f(Sik )  =  (JSi  where,  consid¬ 
ering  the  fact  that  f  is  a  bijection,  it  would  be  true  if  for  any  (in-- .  ,4*)  £  Sik  we  have 
€  U!ii^.  For  /  defined  as  in  (12)  and  any  €  5ifc  we 

have 

m 

fifli  •  •  •  i  4?  •  •  *  j  fm)  =  4  +  7Tj— 1  (4  —  1)  H"  TTfc  — 1  (4  —  1) 

i=2 

m 

=  /(®lj  •  •  •  >4-l)  +  TTfc-lfe  —  1)  +  ^  7Tj_i(4  1)- 

The  ranges  for  f  (iu  •  •  •  ,4-i)  are  [42, .  • . ,  7rfc— i]  and  [0,7r*,  27rfc, . .  .,7rfc+i-~ 

7rjb,7r*+i, . . .  ,7r  —  7 Tfe]  respectively ,  while  the  term  7Tk-i(ik  —  1)  is  constant.  Combining  these 
ranges  we  can  see  that  /(4, . . .  ,4?  ■  •  •  ^ m )  €  . . . , Blukk ] .  □ 

It  follows  from  lemma  4.2  that  we  can  generate  B\k  and  therefore  Sik  in  O(^) 
time. 


5  An  Algorithm  for  the  MAP 

The  algorithm  which  we  will  use  for  solving  the  MAP(ra)  follows  the  general  methodology 
of  Greedy  Randomized  Adaptive  Search  Procedures  (GRASP).  GRASP  is  an  iterative  pro¬ 
cedure  consisting  of  two  phases  at  each  iteration,  a  construction  phase  and  a  local  search 
phase.  In  the  first  phase,  a  solution  is  constructed  in  a  greedy  randomized  fashion,  while  in 
the  second  phase,  a  local  search  is  applied  in  the  neighborhood  of  the  constructed  solution, 
for  possible  further  improvement.  This  process  is  repeated  until  a  stopping  criterion,  such 
as  maximum  number  of  iterations,  is  satisfied.  Each  iteration  can  be  viewed  as  a  procedure 
which  samples  high-quality  points  from  the  solution  space,  and  applying  a  local  search  to 
the  point  chosen  to  obtain  a  local  optimum  For  a  tutorial  and  survey  of  GRASP,  sec  Fco 
and  Resende  [2].  Figure  4  shows  the  pseudo-code  for  a  generic  GRASP. 

GRASP  has  been  implemented  to  solve  Quadratic  Assignment  Problems  (QAP)s  with 
dense  [7,  5]  and  sparse  coefficient  matrices  [11],  while  an  implementation  for  solving  the 
BiQuadratic  Assignment  Problem  can  be  found  in  [8].  A  parallel  implementation  of  a 
GRASP  for  solving  the  QAP  is  described  in  [10].  In  the  sections  that  follow,  the  construction 
phase  and  the  local  search  phase  of  the  algorithm  as  applied  to  the  MAP(m)  are  described. 
It  is  noted  that  the  formulation  used  for  the  development  of  the  algorithm,  is  the  one 
presented  on  §2.2  using  permutations. 
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procedure  GRASP(ListSize , Maxlter , RandomSeed) 

1  InputlnstanceO; 

2  do  k  =  1, ... ,  Maxlter  —* 

3  ConstructSolut ion (ListSize , RandomSeed) ; 

4  LocalSearch(BestSolutionFound); 

5  UpdateSolution(BestSolutionFound) ; 

6  od; 

7  return  BestSolutionFound 
end  GRASP; 


Figure  4:  Pseudo-code  for  a  generic  GRASP 


5.1  Construction  of  a  Solution 

The  objective  function  of  the  MAP  is  nothing  else  but  the  summation  of  m  elements  of 
C  such  that  the  permutation  constraints  are  satisfied.  In  the  construction  phase  we  will 
construct  a  solution  in  a  greedy  randomized  manner.  Let  us  define  the  set  of  already  made 
assignments  upon  making  the  kth  assignment  as 

r k:={ieN  :  f~l(i)  =  (Pi(?)>  ■  •  •  >PmU))  for  somc  j  =  1,  • .  • ,  k  -  i}, 

and  the  set  of  infeasible  assignments  resulting  from  assigning  s  in  the  solution 

M 

As:—{ieN  :  for  /^(s)  =  (ils . . .  ,iM),  *€  |J  S^}. 

j— i 

Our  RCL  will  be  represented  by  a  combination  of  a  doubly  linked  list,  with  a  pointer  array. 
The  array  keyO  is  the  c  array  sorted  in  increasing  order,  and  the  inv()  array  is  just  a 
pointer  array  that  shows  the  position  in  the  keyO  of  any  index  of  c.  The  arrays  prev() 
and  nextO  arc  used  for  the  doubly  linked  list,  that  is  to  move  in  the  keyO  array. 


Figure  5:  Doubly  linked  lists 


Example  5.1  Consider  for  example  Figure  5  where  the  c  array  and  the  invO  ,  prevO  , 
keyO  and  nextO  arrays  are  as  shown .  Assume  that  we  make  the  first  assignment ,  and. 
choose  randomly  among  the  first  [a7rj  smallest  elements  of  keyO,  and  assume  that  we  pick 
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procedure  Construct ionPhase  (a) 

1  T  =  tt  :  r0  =  0; 

2  do  k  =  1, . . .  ,riA/  — ► 

3  t  =random[l,  [aT\]  J 

4  s=startkey; 

5  do  i  =  2, . . . ,  t  — ► 

6  5  =next  (s)  ; 

7  od; 

8  s  =key(s) ; 

9  Tk  =  Tk~i  U  {s}; 

10  GenerateDelta(s,  As)  ; 

11  for  i  e  A5  — ► 

12  T  =  T  —  1 ; 

13  UpdateLists (i, prev, next , key)  ; 

14  As  =  As/{i}; 

15  rof; 

16  od; 

17  return  TnM 

end  ConstructionPhase ; 


Figure  6:  Pseudo-code  for  GRASP  construction  phase 


the  element  key (2)— Then  T\  =  {2}  cmd  assume  that  A2  =  {2,  5},  where  A2  75  generated 
as  described  previously.  We  update  our  RCL  as  follows: 

For  element  2, 
key (inv(2) )=-l , 

prev (next (inv (2) ) ) =prev (inv (2) ) , 
next (prev ( inv (2) ) ) =next ( inv (2) ) , 
and  for  element  5, 
key(inv(5))=-l, 

prev (next ( inv ( 5) ) ) =pre v ( inv (5) ) , 
next (prev ( inv (5) ) ) =next (inv (5) ) . 

Every  time  we  choose  an  element  from  the  key  ( )  array  we  proceed  using  the  prev  ( ) ,  next  ( ) 
arrays. 

The  construction  phase  is  shown  in  pseudo-code  in  Figure  6.  In  line  1  we  initialize  the 
size  or  our  RCL,  and  the  set  T.  In  the  loop  defined  by  the  lines  2-16,  all  the  assignments 
arc  made.  After  generating  a  random  number  in  line  3,  in  lines  4-8  we  move  in  the  doubly 
linked  list  to  retrieve  the  index  s  which  corresponds  to  our  current  assignment.  We  update 
T  in  line  9,  and  in  line  10  we  generate  the  set  of  infeasible  assignments  due  to  s,  the  set  As. 
The  procedure  GenerateDelta(s,  As)  is  done  as  described  in  4.2.  In  lines  11-15  we  update 
the  RCL,  by  excluding  from  the  doubly  linked  list  the  infeasible  elements  contained  in  As. 
The  procedure  UpdateLists(7,prev, next , key)  is  described  earlier. 

5.2  Local  Search 

After  completion  of  the  construction  phase,  we  have  a  feasible  solution  which  is  consisted  of 
m  n\ -permutations.  In  the  local  search  phase,  we  try  to  improve  each  solution  constructed 
in  phase  1  by  searching  its  neighborhood  for  a  better  solution. 

Let  p  be  a  r-permutation  of  n  elements  (r  <  n).  Define  the  difference  between  p  and  a 
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distinct  r-pcrmutation  q,  to  be 

Sr{p,q)-={i  ■  p(i)  +  (13) 

Define  the  distance  between  them  as 

dr(p,q)  ■=  (14) 

The  k-exchange  neighborhood  of  p  is  the  set  of  all  permutations  q  such  that  dr(p,  q)  <  fc, 
Nk(p)  :=  {q  :  dr(p,q)  <  fc,  1  <  A:  <  r}.  (15) 

Remark  5.1  | iVg (p) |  =  Q  +  r(n  -  r). 

For  a  given  P  G  P7n  define  its  2- exchange  neighborhood  as  the  set 

N2{P)  :=  {Q  €  Pm  :  dm(pj,qj)  <  2  for  all  j  =  1, . . .  ,m}.  (16) 

Remark  5.2  |A^(P)|  =  ]ljli  ((”2*)  +  nm(^j  - 

A  local  search  algorithm  based  on  this  neighborhood  definition  would  run  in  exponential 
time,  therefore  we  have  to  employ  a  different  neighborhood  definition  which  is  smaller  in 
size.  Let 


N2{P)  N2{P)  :=  {Q  e  Pm  :  dm(pj,qj)  <  2  for  some  j  =  1, . . .  ,m).  (17) 

Therefore  each  Q  G  N2(P)  differs  from  P  in  only  one  column. 

Remark  5.3  \N2{P)\  —  m(n™)  +  nm  Y?jLi(nj  ~~  nrn)- 

For  every  Q  G  7V2(P)  we  can  associate  a  triplet  ( j ,  fc,  l)  where  j  =  1, . . .  ,m,  k  =  1, . . . ,  nm 
and  l  —  which  indicates  that  pj(k)  =  ^(Z)  and  Pj(Z)  =  qj{k).  If  Z  <  nrn  then 

dm (pj ? Qj )  =  2  while  if  nm  <  l  <  nj  then  dm(pj,  ^)  =  1. 

Next  we  define  the  gain  (either  positive  or  negative)  between  two  solutions  in  a  neigh¬ 
borhood.  Let 

ni 

Z(P)  =  Cjp2(»)-“Pm-l(0pm(0 

for  some  P  €  Prn.  Note  that  Z(P)  is  equivalent  with  the  objective  function  defined  in  (10). 
Then  it  is  immediate  that  for  any  Q  G  N2(P) 

Z(P)  —  Z(Q)  =  C'p1(i)...prn(i)  ~  Oqi  (i)- -9m  (O’ 


and  due  to  the  definition  of  iV2(P)  there  is  only  one  j  such  that  S(pj,qj)  /  0.  Having 
defined  (17)  and  (19),  a  local  search  procedure  for  finding  local  optima  is  straightforward. 
Given  any  P  we  generate  the  gains  Z(P)  -  Z(Q)  for  every  Q  G  iV2(P),  and  we  pick  the  Q 
with  the  largest  positive  gain.  The  procedure  is  repeated  in  N2(Q)  until  no  positive  gains 
can  be  found,  that  is,  our  current  solution  is  the  locally  optimum. 
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6  Problem  Generator 

The  efficiency  of  an  algorithm  for  a  given  combinatorial  optimization  problem  that  belongs 
to  the  NP-Hard  complexity  class,  is  determined  by  several  criteria  including  the  accuracy  of 
the  solution,  the  speed  of  the  algorithm,  and  the  effectiveness  of  the  algorithm  with  respect 
to  different  problem  instances.  Since  it  is  unlikely  (unless  P=NP)  that  a  polynomial  time 
algorithm  could  be  devised  that  would  solve  the  problem  efficiently,  we  can  only  rely  on 
heuristic  procedures  that  will  provide  good  solutions  in  reasonable  computational  time. 

In  evaluating  the  performance  of  heuristic  algorithms  it  is  crucial  to  have  a  set  of  good 
benchmark  problem  instances,  to  perform  the  computational  experiments.  Basically  there 
arc  two  types  of  benchmark  problems  for  a  given  problem,  the  first  involves  instances  which 
come  from  some  real-world  application,  and  the  second  are  instances  which  are  the  outcome 
of  some  problem  generator  that  was  designed  to  construct  random  instances  for  the 
particular  problem.  While  instances  from  a  real-world  application  arc  desirable,  nevertheless 
they  have  the  dissad vantage  that  they  tend  to  present  trends  in  the  data  related  to  the 
specific  application,  thereby  an  algorithm  which  performs  good  at  a  particular  benchmark 
problem  set  might  not  necessarily  perform  the  same  for  some  other  set  coming  from  a 
different  application.  Moreover,  the  availability  of  instances  from  real-world  applications  is 
rather  limited.  Therefore,  the  design  of  good  problem  generators  for  a  particular  problem 
is  desirable.  Various  characteristics  such  as  instance  size,  hardness,  sparsity  etc.  could  also 
be  implemented  in  instances  generated  as  such. 

6.1  A  Random  Generator 

Generally  it  is  hard  to  construct  a  good  problem  generator  for  a  combinatorial  optimization 
problem,  and  not  so  much  research  has  been  devoted  to  the  subject  [12,  3,  1,  6].  One  may 
conjecture  that  the  problem  of  finding  an  algorithm  that  will  generate  problem  instances 
that  are  NP-hard  for  a  given  problem  class  is  itself  an  NP-hard  problem.  The  characteristics 
that  a  good  problem  generator  should  have  can  be  the  following: 

1.  problem  instance  characteristics  such  as  size,  sparsity  etc.  should  be  a  parameter  of 
the  generator  algorithm. 

2.  the  problem  instances  that  it  generates  should  not  be  easy  to  solve  by  any  heuristic 
or  exact  method. 

3.  the  hardness  of  the  instances  generated  should  be  a  parameter  of  the  generator  algo¬ 
rithm. 

The  first  item  on  the  list  above  is  easy  to  implement.  However  items  2  and  3  are  not  trivial. 

Here  we  propose  an  initial  problem  generator  that  generates  random  instances  of  the 
MAP.  The  algorithm  is  simple,  and  it  proceeds  by  first  generating  a  random  multidimen¬ 
sional  array  C,  where  the  costs  are  random  numbers  distributed  uniformly  in  the  interval 
[2,  M]  where  M  is  a  parameter  of  the  generator.  Then  all  the  diagonal  elements  of  the  array 
arc  set  to  1,  thereby  making  the  optimal  solution  to  be  n\  (i.e.  the  smallest  range),  and  the 
optimal  permutations  are  the  identity  permutations. 

7  Computational  Results 

In  this  section,  computational  results  that  illustrate  the  performance  of  the  algorithm  for 
solving  the  multidimensional  assignment  problem  are  presented.  The  algorithm  imple¬ 
mented  follows  the  description  given  in  §5,  while  the  test  problems  were  generated  using 
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algorithm  ConstructRandomlnst ance  (T,  m,  ni ,  , . . 

■  j  ) 

1 

*:=n£i»i 

2 

do  i  =  1,2, . . .  ,7r  — * 

3 

(U  5  ^2?  •  •  •  5  Vn)  ~  /  (0  j 

4 

C(«„i2 . im)  =  random[2,T]; 

5 

od; 

6 

do  i  =  1,2, . . .  ,ni  — ► 

7 

=  i; 

8 

od; 

9 

return  C\ 

end  ConstructSparse; 

Figure  7:  MAP  random  instance  generator 


the  problem  generator  described  in  §6.  Overall  twenty  one  problems  were  generated  of  var¬ 
ious  sizes  and  dimension.  For  each  problem  the  GRASP  algorithm  was  allowed  to  run  for  a 
maximum  of  100000  iterations  or  until  the  optimum  solution  was  found.  The  computational 
runs  were  performed  on  a  Intel  Pentium  II  400  MHz  computer,  with  128MBytes  of  memory, 
running  the  Linux  operating  system.  For  each  problem  instance  the  number  of  iterations 
and  the  CPU  time  until  the  algorithm  found  the  optimum  solution  was  recorded.  If  the 
optimum  was  not  found,  the  aforementioned  values  correspond  to  the  best  solution  found 
by  the  algorithm. 

From  the  computational  results  in  Table  1  we  can  see  that  the  algorithm  finds  the 
optimum  solution  in  all  but  three  instances.  The  instances  for  which  the  optimum  is  not 
found  within  the  maximum  number  of  iterations,  arc  amongst  the  largest  ones  in  terms  of 
number  of  elements. 


8  Conclusions  and  Future  Research 

The  primary  subject  of  study  in  this  research  is  the  multidimensional  assignment  prob¬ 
lem  (MAP)  as  it  arises  from  its  application  in  a  multiple-sensor  multiple-target  tracking 
framework.  The  MAP  is  an  NP-Hard  combinatorial  optimization  problem,  which  has  been 
studied  before  in  the  past.  It  is  a  generalized  version  of  the  well  known  linear  assignment 
problem,  and  it  seems  to  enjoy  various  applications  from  different  fields. 

In  this  research  we  provide  new  mathematical  formulations  of  the  problem,  all  of  them 
equivalent  to  each  other,  and  to  the  classical  integer  programming  formulation.  Using  the 
permutation  formulation  of  the  problem,  which  has  an  inherent  combinatorial  nature,  a 
purely  discrete  heuristic  algorithm  for  providing  close  to  optima  solutions  to  the  problem  is 
presented.  In  order  for  the  algorithm  to  perform  satisfactory  fast,  efficient  data  structures 
to  handle  the  massive  data  associated  with  the  problem  arc  developed,  and  integrated  to  the 
algorithm.  The  result  is  a  heuristic  algorithm  that  can  solve  MAP’s  of  varying  dimension. 
Moreover,  using  the  same  permutation  formulation  of  the  problem,  a  problem  generator  that 
generates  MAP  instances  of  known  optimal  solution  is  developed  and  implemented.  The 
instances  generated  arc  purely  for  testing  purposes.  All  algorithms  have  been  implemented 
in  FORTRAN,  and  copies  of  working  codes  are  provided. 

An  interesting  aspect  that  emerged  in  the  process  of  this  research,  is  the  formulation  of 
the  problem  as  a  matroid  intersection  problem.  Using  results  from  the  rich  matroid  theory, 
this  new  approach  might  lead  to  new  and  power  full  solution  methods  for  the  problem. 
Specifically  the  following  open  questions  emerged: 
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GRASP  Solution 


name 

dimension 

size 

optimum 

iterations 

CPU  time 

solution 

6x4 

6 

4096 

4 

364 

Omls 

4 

4x10 

4 

10000 

10 

166 

0m2s 

10 

3x22 

3 

10648 

22 

129 

0m05s 

22 

6x5 

6 

15625 

5 

2228 

0m35s 

5 

5x7 

5 

16807 

7 

7234 

2ml5s 

7 

7x4 

7 

16384 

4 

2204 

0m36s 

4 

5x8 

5 

32768 

8 

1805 

lm23s 

8 

4x14 

4 

38416 

14 

288 

0ml  9s 

14 

6x6 

6 

46656 

6 

2627 

3m  16s 

6 

4x15 

4 

50625 

15 

5078 

7m51s 

15 

3x38 

3 

54872 

38 

176 

0m53s 

38 

5x9 

5 

59049 

9 

9228 

16m04s 

9 

3x40 

3 

64000 

40 

149 

0m37s 

40 

8x4 

8 

65536 

4 

5470 

llm30s 

4 

7x5 

7 

78125 

5 

2228 

5m37s 

5 

3x46 

3 

97336 

46 

1772 

15m55s 

46 

5x10 

5 

100000 

10 

1477 

4m33s 

10 

6x7 

6 

117649 

7 

7651 

30m40s 

15* 

4x20 

4 

160000 

20 

4876 

26ml0s 

20 

7x6 

7 

279936 

6 

3226 

10m22s 

12* 

8x5 

8 

390625 

5 

662 

15ml0s 

13* 

Table  1:  Summary  of  computational  runs  (*  indicates  not  optimal) 


•  How  docs  the  dual  problem  structure  looks  like,  and  if  there  is  a  possibility  to  develop 
a  primal-dual  algorithm  for  solving  the  problem.  Such  primal-dual  algorithms  on 
matroid  intersections  that  involve  two  matroids  are  optimal,  but  in  the  case  of  MAP 
where  we  have  the  intersection  of  more  than  two  matroids,  an  approximation  algorithm 
could  be  developed. 

•  Is  there  an  efficient  way  to  compute  or  approximate  (bound)  the  value  of  p(E ,  !F)1  If 
such  an  algorithm  is  constructed,  then  we  immediately  have  an  approximation  algo¬ 
rithm  for  the  minimazation  version  of  the  MAP(m).  Note  that  for  the  maximization 
version  of  the  problem,  a  similar  approach  gives  an  e- approximation  algorithm  with 
e  =  1/m. 

m 

•  Develop  a  basis  oracle  for  the  independent  system  (E,  p|  Tv3 ) •  That  is,  an  algorithm 

i=1 

that  given  some  ICE  will  determine  whether  or  not  I  contains  a  maximal  independent 
set. 

•  Develop  a  problem  generator  using  the  matroid  intersection  formulation. 

9  FORTRAN  Codes 

9.1  Problem  Generator 

c  - 

c  Program  to  generate  random  input  data  for  n-dimensional 

c  assignment  problem.  The  elements  of  the  cost  matrix  C  are 
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c  between  1-100.  The  problem  instances  have  known  optimal 
c  solutions,  the  identity  permutations,  with  optimal  solution  value 
c  equal  to  Nmin  :=  the  minimum  of  Nl,...,Ndim. 
c 

c  Author:  Pitsoulis  Leonidas 

c  - 

c  Inputs  from  user: 

c 

c  dim  -  the  dimension 

c  Nl, .  .  . ,Ndim  -  the  range  for  each  index  1. .dim 

c 

integer  nmax,dimmax,npower 
parameter  (nmax=100 ,dimmax=10) 
parameter  (npower=10000000) 

integer  c(npower),  nindex(dimmax) ,  indx(dimmax) 
integer  dim,  n 
integer  seed, high, i,  j,k 
real  randp,xrand 

c  - 

seed=2700001 

high=100 

print* , ’Enter  dimension  dim:’ 
read* ,  dim 

print*, *  Enter  Nl, . . . ,Ndim  (in  ascending  order) : * 

read*,  (nindex(i),  i=l,dim) 

n=l 

do  10  i=l,dim 
n=n*nindex(i) 

10  continue 

do  20  i-l,n 

xrand  =  randp(seed) 
nselct  =  1  +  seed/ (2147483647/high) 
c(i)  =  nselct  +  1 
20  continue 

do  30  i=l ,nindex(l) 
do  40  k=l ,dim 
indx(k)  =  i 
40  continue 

call  indxi (n , dim , indx , j , nindex) 
c(j)  =  1 
30  continue 

do  50  i=l,n 

print  *,  c(i) 

50  continue 


end 

real  function  randp(ix) 

c  - 

c  randp:  Portable  pseudo-random  number  generator, 
c  Reference:  L.  Schrage,  "A  More  Portable  Fortran 
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c  Random  Number  Generator",  ACM  Transactions  on 

c  Mathematical  Software,  Vol.  2,  No.  2,  (June,  1979). 

c  - 

integer  a,p, ix,bl5 ,bl6 ,xhi ,xalo , leftlo ,fhi ,k 
data  a/16807/, M5/32768/ ,bl6/65536/ , p/2147483647/ 

c  - 

xhi  =  ix/bl6 

xalo  =  (ix  -  xhi*bl6)*a 

leftlo  =  xalo/bl6 

fhi  =  xhi*a  +  leftlo 

k  =  fhi/bl5 

ix  =  (((xalo  -  leftlo*bl6)  -  p)  +  (fhi  -  k*bl5)*bl6)  +  k 
if  (ix  .It.  0)  ix  =  ix  +  p 

randp  =  float (ix)*4.656612875e-10 

c  - 

return 

c  end  of  randp 
end 


9.2  Data  Structures 

9.2.1  The  Cost  Array  Functions 

subroutine  iindx (n , dim , i , nindex , indx) 


Given  the  dimension  dim  of  the  multi-dim  array  C,  which  is 
represented  as  1-dim  array  c(l)...c(n),  then  given  the  index  i 
of  some  element  c(i)  return  the  dim-indices  (dim=2 ,3, . . etc . )  that 
correspond,  in  the  array  indx() . 

integer  nindex (dim) , indx (dim) 
integer  n,dim,i,j,l 
real  p 
l=i 

p=real (n) 
do  10  j=dim,l,-l 

p=p/ (real (nindex ( j ) ) ) 
indx(j)=int ((real(l-l) )/p)  +  1 
1=1-  (indx(j)  -  1  )*int(p) 
continue 
return 
end 

subroutine  indxi (n, dim, indx, j , nindex) 

c  - 

c  Given  the  dim-indices  (dim=2 ,3, . . .etc . )  in  the  array  indx() , 
c  return  the  index  j 

c  - 

integer  nindex (dim) , indx (dim) 
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integer  n,dim,j,k 
real  p 
p=real (n) 
j=0 

do  10  k=dim,l,“l 

p=p/ (real (nindex (k) ) ) 
j=j+p*(indx(k)-l) 

10  continue 

j*j+l 

return 

end 


9.2.2  The  Queues 

subroutine  insrtq(n,v,iv, sizeq, q,iq) 

c  - 

c  insrtq:  Insert  an  element  (v,iv)  into  a  queue  (q,iq). 

c  - 

c  Passed  scalars: 
c 

c  n  -  size  of  c 

c  v  heap  element  (value) 

c  iv  heap  element  (index) 

c  sizeq  -  size  of  heap 

c 

c  - 

integer  n, sizeq, iv,v 

c  - 

c  Passed  arrays : 

c 

c  q  heap  (value) 

c  iq  heap  (index) 

c 

c  - 

integer  iq(n),q(n) 

c  - 

c  Local  scalars: 

c 

c  sq  temporary  size  of  heap 

c  tsz  -  temporary  variable  (sq/2) 

c  w  -  temporary  variable 

c  iw  temporary  variable 

c 

c  - 

integer  sq,tsz,w,iw 

c  - 

c  Insert  element  into  heap. 

c  - 

sizeq=sizeq+l 
q(sizeq) =v 
iq(sizeq)=iv 

c  - 

c  Update  heap  to  proper  order. 
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sq=sizeq 

w=q(sq) 

iv=iq(sq) 

10  tsz=sq/2 

if  (tsz.ne.O)  then 

if  (q(tsz) .  ge .  v)  then 
q(sq)=q(sq/2) 
iq(sq)=iq(sq/2) 
sq=sq/2 
goto  10 
end  if 
endif 
q(sq)=w 
iq(sq)-iw 
return 
end 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  removq(n , v , iv , sizeq , q , iq) 


removq:  Remove  smallest  element  (v,iv)  from  a  priority- 
queue  (q,  iq)  . 


Passed  scalars: 
n  -  size  total 

v  -  smallest  element  in  heap  (value) 
iv  -  smallest  element  in  heap  (index) 
sizeq  -  size  of  heap 


integer  n,v,iv, sizeq 


Passed  arrays : 

q  -  heap  (value) 

iq  -  heap  (index) 


integer  iq(n),q(n) 


Local 

scalars : 

vtmp 

-  tmp  smallest 

element 

in  heap 

(value) 

ivtmp 

-  tmp  smallest 

element 

in  heap 

(index) 

k 

-  heap  counter 

J 

-  heap  counter 

(2*k) 

szqd2 

-  sizeq/2 

integer  vtmp , ivtmp , k ( j , szqd2 


c 


Remove  element  from  heap. 
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v=q(l) 

iv=iq(l) 

q(l)=q(sizeq) 

iq(l)=iq(sizeq) 

sizeq=sizeq-l 


c  Update  heap  to  proper  order. 


k=l 

vtmp=q(k) 

ivtmp=iq(k) 

szqd2=sizeq/2 

10  if  (k  .gt.  szqd2)  goto  20 

j“k+k 

if  (j  .It.  sizeq)  then 
if  (q(j)  .gt.  q(j+l) )  j-j+i 
endif 

if  (vtmp  .le.  q(j))  goto  20 

q(k)=q(j) 

iq(k)=iq(j) 

k-j 

goto  10 
20  continue 

q(k)=vtmp 
iq(k)=ivtmp 
return 
end 

9.3  Local  Search 

subroutine  local (n , dim , nmax , dimmax , nsmall ,  c , pm , nindex , ob j v , 

+  indexi, indexj ) 

c  - 

c  Perforin  local  search  in  the  solution  pm(.,.)  with  objv  as  the 
c  objective  function  value. 

c  - 

c  - 

c  Passed  scalars  and  arrays : 


integer  dim , dimmax , nmax , nsmall , ob  j  v , n 

integer  c (n) , nindex (dim) , pm (dimmax , nmax) , indexi (dim) , index j (dim) 


c  Local  scalars  and  arrays: 


integer  gain , dif f , i , j , k , 1 , io , j  o , in , j  n , kk , temp , ic , j  c , kc 
real  p 
10  gain=0 

do  20  k=l,dim 

do  30  i-1,  nsmall 

do  40  j«i+l,  nindex (k) 
do  50  1=1, dim 

indexi (l)=pm(l,i) 
indexj (l)=pm(l, j) 

50  continue 
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temp=indexi(k) 

indexi (k) =indexj (k) 

index j (k)-temp 

p=real(n) 

io=0 

jo=0 

in=0 

jn=0 

do  60  kk=dim,l,-l 

p-p/ (real (nindex (kk) ) ) 
io=io+p*(pm(kk, i)-l) 
j  o= j  o+p* (pm (kk , j ) - 1 ) 
in=in+p* (indexi (kk) -1) 
jn=jn+p*( index j (kk)-l) 

60  continue 

io=io+l 
jo=jo+l 
in=in+l 
jn=jn+l 

if  (j . le .nsmall)  then 

diff=  c(io)+c(jo)-c (in)-c(jn) 
else 

diff=  c(io)-c(in) 
endif 

if  (diff .gt .gain)  then 
gain=diff 
ic=i 

jc=j 

kc=k 

endif 

40  continue 

30  continue 

20  continue 

if  (gain.gt.0)  then 
objv=obj v-gain 
temp=pm(kc,ic) 
pm (kc , ic ) =pm (kc , j  c ) 
pm(kc, jc)=temp 
goto  10 
endif 
return 
end 

9.4  Construction  of  a  Greedy  Solution 

subroutine  construct (n , dim , nmax , dimmax , nsmall , c , pm , inv , prev , next , 
+  key , seed , alpha , nindex , indx , batch , numbatch , period, 

+  objv) 


c  Construct  a  greedy  randomized  feasible  solution 


integer  n,dim, nmax, nsmall ,i , j ,k,l,m, objv, dimmax 
integer  c (n) , inv (n) , prev (n) , next (n) , key (n) 
integer  nindex (dim) , indx (dim) , pm (dimmax , nmax) 
integer  batch (dim) , numbatch (dim) , period (dim) 
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integer  high , rclsiz , seed , high , nselct , start , chosen , index , exi , exj 

real*4  r an dp, xr and, alpha 

rclsiz=n 

start=l 

objv=0 

do  10  i=l,nsmall 

c  - 

c  Select  element,  at  random,  from  the  best  (rclsiz) *alpha  cost 
c  elements . 

c  - 

xrand=randp (seed) 
high=alpha* (rclsiz) 
if  (high.lt. 1)  then 
high=l 
endif 

nselct=l+seed/ (2147483647/high) 

c  - 

c  The  variable  chosen  will  be  the  nselct-th  smallest  element  chosen 
c  from  the  double  linked  list,  i.e.  key (chosen) =index  in  c()  of 
c  element  being  chosen. 

c  - 

chosen=start 
do  20  j =2, nselct 

chosen=next (chosen) 

20  continue 

index=key (chosen) 

c  - 

c  Exclude  from  the  doubly-linked  list  all  the  elements  that  are  not 

c  feasible  anymore  due  to  the  insertion  of  the  c (index)  element 

c  in  the  constructed  solution. 

c  - 

call  iindx (n , dim , index , nindex , indx) 
do  25  1=1,  dim 

do  26  m=i,nindex(l) 

if  (pm(l,m) .ne.indx(l))  then 
goto  26 
endif 

temp=pm(l,i) 
pm(l,i)=pm(l,m) 
pm(l,m)=temp 
goto  25 

26  continue 

25  continue 

objv=obj v+c (index) 
do  30  j=l,dim 
exi=(indx(j)-l)*batch(j) 
do  40  k=l,  numbatch(j) 
do  50  m=l,  batch(j) 
exi=exi+l 

c  - 

c  exclude  the  element  exi  from  the  doubly  linked  list 

c  - 

exj=inv(exi) 

if  (key(exj) . eq. -1)  then 
goto  50 
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endif 

key(exj)«-l 
rclsiz=rclsiz-l 
if  (prev(exj) .eq.-l)  then 
start =next (exj) 
prev (next (exj ) ) =-l 
goto  50 
endif 

if  (next (exj) .eq.-l)  then 
next (prev (exj ) ) =-l 
goto  50 
endif 

next (prev (exj ) ) =next (exj ) 
prev(next (exj ) ) =prev (exj ) 

50  continue 

exi= (indx ( j ) -1 ) *batch( j ) +k*period( j ) 

40  continue 

30  continue 

10  continue 

return 
end 

9.5  Main  Driver  of  the  Code 

C  - 

c  Program  driver  for  GRASP  for  the  Multidimensional  Assignment 

c  Problem  code 

c  Author:  Pitsoulis  Leonidas 

c  - 

c  This  file  includes  the  following  Fortran  subroutines  and 
c  functions: 

c 

c  gmap 

c  srtcst 

c  stage 1 

c  stage2 

c  local 

c  insrtq 

c  removq 

c  randp 

c 

c  - 

c  The  input  data  consists  of  the  following  variables  and  arrays: 
c 

c  dim  -  the  dimension  of  the  array 

c  nindexO  -  the  array  with  the  index  ranges  for  each  dimension 

c  nindex(l),  ...  ,  nindex(dim) 

c  c()  -  the  array  elements  c (1) , . . . , c (nindex(l) * . . . *nindex(dim) ) 

c 

c  - 

c  The  array  c  is  represented  as  one-dimensional  array  c(l)...c(n) 
c  where  the  variable  n  =  nindex(l) * . . . *nindex(dim) . 

c  Given  and  index  i  in  {l...n}  we  can  retrieve  the  individual 

c  indices  represented  in  the  variable  array  indx(dim), 
c  indx(l),  indx(2),  ...  ,  indx(dim)  using  the  subroutine 


-  control  subroutine  for  GRASP  for  MAP  algorithm 

-  sorts  cost 

-  stage  1  of  GRASP  construction  phase 

-  stage  2  of  GRASP  construction  phase 

-  2-exchange  local  search  for  QAP 

-  insert  element  into  heap  for  sorting 

-  remove  element  from  heap 

-  random  number  generator  function 
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c  iindx ( — ).  Given  indx ( 1) , . . . , indx(dim)  we  can  retrieve  the 

c  index  i  in  {l...n}  using  the  subroutine  indxi( — ). 

c  - 

integer  nmax,  dimmax,  npower 
parameter  (nmax=22 ,dimmax=5) 
parameter  (npower=15000) 
integer  c(npower) 

integer  srtc (npower ) , srtic (npower) , csrt (npower) , icsrt (npower) 
integer  inv(npower) ,prev(npower) ,next (npower) ,key (npower) 
integer  nindex (dimmax) , indx (dimmax) , indexi (dimmax) , index j (dimmax) 
integer  pm (dimmax, nmax) 

integer  batch(dimmax) ,numbatch (dimmax) .period (dimmax) 

integer  dim, n, nsmall 

integer  seed, objv, iter , i ,  j , 1 .totobjv 

real *4  alpha 

real*4  cput(2) ,bcpu,etime,totcpu 

c  - 

C  Read  Problem  input  data  and  parameters 

c  - 

call  r eadp (nmax , dimmax , npower , dim , c , nindex , n , nsmall ) 

c  - 

C  Sort  the  array  c()  into  csrt()  to  be  used  and  the  corresponding 
C  indices  of  csrt()  into  icsrt () 

c  - 

call  srtcst (n , dim, c , srtc , srtic , csrt , icsrt) 

c  - 

c  Construct  the  array  inv() 

c  - 

call  const inv(n, dim, inv, icsrt) 

call  initial (n , dim , dimmax , pm , nmax , ob j v , batch , numbat ch , period , 

+  nindex) 

seed=270001 
alpha=0 . 1 
iter=100 
objv=0 

totobjv=9999999 

c  - 

c  Main  GRASP  iterations. 

c  - 

bcpu=etime (cput) 
do  10  i=l,iter 

call  build(n,prev,next , icsrt , key) 

call  construct (n , dim , nmax , dimmax , nsmall , c , pm , inv , prev , next , key , 
■f  seed ,  alpha ,  nindex ,  indx ,  bat  ch ,  numbat  ch ,  period ,  objv) 

call  local (n , dim , nmax , dimmax , nsmall , c , pm , nindex , ob j v , indexi , 

+  index j) 

c  call  chkcst (n, dim, nmax, dimmax, pm, c, nindex, costos, nsmall) 

if  (objv.lt .totobjv)  then 

print  *,  1 #######  at  iteration  :  ’  .i,*  ######### } 
do  20  1=1, dim 

print  *,  (pm(l , j ) ,  j=l , nindex (1) ) 

20  continue 

print  *,  * - Cost  is - >  objv 

c  print  *,  * - Acutal  Cost  is - >  * ,  costos 


10  BIBLIOGRAPHY 


29 


totobjv^objv 
end  if 
10  continue 

totcpu=etime (cput) -bcpu 

print  *,  * - time - >  * ,  totcpu 

stop 

end 


10  Bibliography 
References 

[1]  K.  Bassam,  P.M.  Pardalos,  and  D.Z.  Du.  A  test  problem  generator  for  the  stcincr 
problem  in  graphs.  ACM  Transactions  on  Mathematical  Software ,  19(4) :509— 522,  1993. 

[2]  T.A.  Fco  and  M.G.C.  Resende.  Greedy  randomized  adaptive  search  procedures.  Journal 
of  Global  Optimization,  6:109-133,  1995. 

[3]  J.  Hasselbcrg,  P.M.  Pardalos,  and  G.  Vairaktarakis.  Test  case  generators  and  com¬ 
putational  results  for  the  maximum  clique  problem.  Journal  of  Global  Optimization , 
3:463-482,  1993. 

[4]  B.  Korte  and  C.L.  Monma.  Some  remarks  on  a  classification  of  oracle- type  algorithms. 
In  Numerische  Methoden  bei  graphentheoretischen  und  kombinatorischen  Problemen, 
Band  2 ,  pages  195-215.  Birkhauser,  Basel,  1979. 

[5]  Y.  Li,  P.  M.  Pardalos,  and  M.  G.  C.  Resende.  Algorithm  754:  FORTRAN  subroutines 
for  approximate  solution  of  dense  quadratic  assignment  problems  using  GRASP.  ACM 
Transactions  on  Mathematical  Software ,  22:104-118,  1996. 

[6]  Y.  Li  and  P.M.  Pardalos.  Generating  quadratic  assignment  problems  with  known  opti¬ 
mal  permutations.  Computational  Optimization  and  Applications ,  1(2):  163-184,  1992. 

[7]  Y.  Li,  P.M.  Pardalos,  and  M.G.C.  Resende.  A  greedy  randomized  adaptive  search 
procedure  for  the  quadratic  assignment  problem.  In  P.M.  Pardalos  and  H.  Wolkowicz, 
editors,  Quadratic  assignment  and  related  problems ,  volume  16  of  DIM  ACS  Series  on 
Discrete  Mathematics  and  Theoretical  Computer  Science ,  pages  237-261.  American 
Mathematical  Society,  1994. 

[8]  T.  Mavridou,  P.  M.  Pardalos,  L.  S.  Pitsoulis,  and  M.  G.  C.  Resende.  A  GRASP 
for  the  biquadratic  assignment  problem.  European  Journal  of  Operations  Research , 
105(3):613-621,  1998. 

[9]  R.  A.  Murphey,  P.  M.  Pardalos,  and  L.  Pitsoulis.  A  greedy  randomized  adaptive  search 
procedure  for  the  multitarget  multisensor  tracking  problem.  In  Network  Design:  Con¬ 
nectivity  and  Facility  Location,  volume  40  of  DIM  ACS  Series  on  Discrete  Mathematics 
and  Theoretical  Computer  Science,  pages  277-302.  American  Mathematical  Society, 
1997. 

[10]  P.  M.  Pardalos,  L.  S.  Pitsoulis,  and  M.  G.  C.  Resende.  A  parallel  GRASP  implementa¬ 
tion  for  solving  the  quadratic  assignment  problem.  In  A.  Ferreira  and  Jose  D.P.  Rolim, 
editors,  Parallel  Algorithms  for  Irregular  Problems:  State  of  the  Art,  pages  115-133. 
Kluwer  Academic  Publishers,  1995. 


REFERENCES 


30 


[11]  P.  M.  Pardalos,  L.  S.  Pitsoulis,  and  M.  G.  C.  Rcsendc.  Algorithm  769:  FORTRAN 
subroutines  for  approximate  solution  of  sparse  quadratic  assignment  problems.  ACM 
Tmnsactions  on  Mathematical  Software ,  23:196-208,  1997. 

[12]  P.M.  Pardalos.  Generation  of  large-scale  quadratic  programs  for  use  as  global  opti¬ 
mization  test  problems.  ACM  Transactions  on  Mathematical  Software ,  13(2):133-137, 
1987. 

[13]  A.  B.  Poore  and  N.  Rijavcc.  Partitioning  multiple  data  sets:  multidimensional  assign¬ 
ments  and  lagrangian  relaxation.  In  P.M  Pardalos  and  H.  Wolkowicz,  editors,  Quadratic 
assignment  and  related  problems ,  volume  16  of  DIM  ACS  Series  on  Discrete  Mathemat¬ 
ics  and  Theoretical  Computer  Science ,  pages  317-342.  American  Mathematical  Society, 
1994. 

[14]  D.  J.  A.  Welsh.  Matroid  Theory .  Academic  Press,  1976. 


